Chapter 8 HMSC analysis

8.1 Load data

load("data/data.Rdata")
load("hmsc/fit_model1_250_10.Rdata")

8.2 Variance partitioning

# Compute variance partitioning
varpart=computeVariancePartitioning(m)

varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
   group_by(variable) %>%
   summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
   tt()
tinytable_y3bzzqyttyssdfioa85l
variable mean sd
Random: location 37.900015 25.317903
logseqdepth 56.110626 25.796874
sex 4.937460 5.612719
origin 1.051899 1.282563
# Basal tree
varpart_tree <- genome_tree

#Varpart table
varpart_table <- varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(genome=factor(genome, levels=rev(varpart_tree$tip.label))) %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location"))))

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
varpart_tree <- varpart_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************
# Add phylum colors next to the tree tips
varpart_tree <- gheatmap(varpart_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
   scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
varpart_tree <- varpart_tree + new_scale_fill()

# Add variance stacked barplot
vertical_tree <-  varpart_tree +
       scale_fill_manual(values=c("#506a96","#cccccc","#be3e2b","#f6de6c"))+
        geom_fruit(
             data=varpart_table,
             geom=geom_bar,
             mapping = aes(x=value, y=genome, fill=variable, group=variable),
             pwidth = 2,
             offset = 0.05,
             width= 1,
             orientation="y",
             stat="identity")+
      labs(fill="Variable")

vertical_tree

8.3 Posterior estimates

# Select desired support threshold
support=0.9
negsupport=1-support

# Basal tree
postestimates_tree <- genome_tree

# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
    mutate(value = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    #select(genome,sp_vulgaris,area_semi,area_urban,sp_vulgarisxarea_semi,sp_vulgarisxarea_urban,season_spring,season_winter,sp_vulgarisxseason_spring,sp_vulgarisxseason_winter) %>%
    column_to_rownames(var="genome")

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)


colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
postestimates_tree <- postestimates_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************
#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
      scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

# Add posterior significant heatmap

postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.5, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
        scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
        labs(fill="Trend")

postestimates_tree +
        vexpand(.25, 1) # expand top 

8.3.1 Positively associated genomes

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Positive") %>%
    arrange(-value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum, class, family)%>%
    tt()
tinytable_l1xwr6tidyexw4ryx4pp
genome phylum class order family species value
bin_CaboVerde.50 Actinomycetota Actinomycetia Actinomycetales Actinomycetaceae NA 0.953
bin_Aruba.25 Actinomycetota Actinomycetia Actinomycetales Bifidobacteriaceae Bifidobacterium pseudocatenulatum 0.995
bin_Aruba.2 Actinomycetota Actinomycetia Actinomycetales Bifidobacteriaceae Bifidobacterium adolescentis 0.994
bin_Aruba.6 Actinomycetota Actinomycetia Actinomycetales Bifidobacteriaceae Bifidobacterium longum 0.985
bin_Denmark.14 Actinomycetota Actinomycetia Actinomycetales Bifidobacteriaceae Bifidobacterium gallinarum 0.950
bin_Aruba.47 Actinomycetota Actinomycetia Mycobacteriales Mycobacteriaceae NA 0.964
bin_Aruba.57 Actinomycetota Actinomycetia Mycobacteriales Mycobacteriaceae Corynebacterium pyruviciproducens 0.931
bin_Aruba.51 Actinomycetota Coriobacteriia Coriobacteriales Atopobiaceae Parolsenella sp900544655 0.969
bin_Denmark.44 Actinomycetota Coriobacteriia Coriobacteriales Atopobiaceae UBA7748 sp900314535 0.951
bin_Aruba.14 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella stercoris 1.000
bin_Denmark.4 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella stercoris 1.000
bin_Aruba.21 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae NA 0.999
bin_Brazil.61 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella tanakaei 0.998
bin_Aruba.39 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella sp902362275 0.992
bin_Spain.84 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella sp900555555 0.992
bin_Brazil.151 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella intestinalis 0.991
bin_Denmark.127 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella sp900548365 0.991
bin_Denmark.17 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella bouchesdurhonensis 0.983
bin_Brazil.81 Bacillota Bacilli RFN20 CAG-826 NA 0.907
bin_Brazil.109 Bacillota_A Clostridia Oscillospirales Butyricicoccaceae AM07-15 sp003477405 0.936
bin_Aruba.28 Bacillota_A Clostridia Oscillospirales Oscillospiraceae NA 0.995
bin_Malaysia.9 Bacillota_A Clostridia Oscillospirales Oscillospiraceae Dysosmobacter welbionis 0.986
bin_Malaysia.116 Bacillota_A Clostridia Oscillospirales Oscillospiraceae Flavonifractor plautii 0.981
bin_Malaysia.34 Bacillota_A Clostridia Oscillospirales Oscillospiraceae Lawsonibacter sp000177015 0.954
bin_Aruba.54 Bacillota_A Clostridia Tissierellales Peptoniphilaceae NA 0.959
bin_CaboVerde.35 Bacillota_A Clostridia Tissierellales Peptoniphilaceae NA 0.930
bin_Aruba.36 Bacillota_A Clostridia Tissierellales Peptoniphilaceae Peptoniphilus_C sp902363535 0.913
bin_Aruba.52 Bacillota_A Clostridia Tissierellales Peptoniphilaceae NA 0.906
bin_Malaysia.26 Bacillota_A Clostridia Oscillospirales Ruminococcaceae NA 0.976
bin_Aruba.29 Bacillota_A Clostridia Oscillospirales Ruminococcaceae Gemmiger variabilis_C 0.941
bin_Malaysia.103 Bacillota_C Negativicutes Acidaminococcales Acidaminococcaceae Acidaminococcus intestini 0.942
bin_Malaysia.8 Bacillota_C Negativicutes Veillonellales Dialisteraceae Dialister sp900543165 0.981
bin_Denmark.60 Bacillota_C Negativicutes Veillonellales Dialisteraceae Allisonella histaminiformans 0.938
bin_Aruba.64 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Megasphaera sp000417505 0.970
bin_Brazil.58 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Megasphaera elsdenii 0.969
bin_CaboVerde.38 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Megasphaera sp000417505 0.961
bin_Malaysia.81 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Megasphaera stantonii 0.943
bin_Malaysia.96 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Caecibacter sp003467125 0.919
bin_Brazil.163 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella lascolaii 0.992
bin_CaboVerde.18 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella copri 0.981
bin_Aruba.10 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella sp900544825 0.978
bin_Brazil.5 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella sp900540415 0.929
bin_Spain.21 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola coprophilus 0.929
bin_Malaysia.18 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotellamassilia sp000437675 0.924
bin_Malaysia.117 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotellamassilia sp900541335 0.919
bin_Malaysia.77 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola sp900542985 0.913
bin_Brazil.75 Campylobacterota Campylobacteria Campylobacterales Helicobacteraceae NA 0.964
bin_Malaysia.61 Campylobacterota Campylobacteria Campylobacterales Helicobacteraceae Helicobacter_C labetoulli 0.964
bin_Brazil.128 Campylobacterota Campylobacteria Campylobacterales Helicobacteraceae NA 0.963
bin_Brazil.70 Desulfobacterota Desulfovibrionia Desulfovibrionales Desulfovibrionaceae Mailhella sp003150275 0.944
bin_Brazil.146 Pseudomonadota Alphaproteobacteria RF32 CAG-239 CAG-495 sp001917125 0.969
bin_Brazil.9 Pseudomonadota Alphaproteobacteria RF32 CAG-239 CAG-495 sp000436375 0.919
bin_Aruba.15 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae NA 0.963
bin_Brazil.51 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae Sutterella wadsworthensis_A 0.963
bin_Spain.43 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae CAG-521 sp000437635 0.959
bin_Brazil.64 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae NA 0.948
bin_Brazil.92 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae CAG-521 sp000437635 0.948
bin_Brazil.15 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae NA 0.930
bin_CaboVerde.33 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Anaerobiospirillum sp900543125 0.995
bin_Brazil.82 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Anaerobiospirillum succiniciproducens 0.990
bin_CaboVerde.55 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Anaerobiospirillum_A thomasii 0.939

8.3.2 Negatively associated genomes

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Negative") %>%
    arrange(value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum,class,family)%>%
    tt()
tinytable_bbzfzysx61kn9u9o37iv
genome phylum class order family species value
bin_Denmark.27 Bacillota Bacilli Lactobacillales Enterococcaceae Enterococcus_B hirae 0.000
bin_Brazil.12 Bacillota Bacilli Lactobacillales Enterococcaceae Enterococcus faecalis 0.001
bin_Brazil.170 Bacillota Bacilli Lactobacillales Enterococcaceae NA 0.014
bin_CaboVerde.24 Bacillota Bacilli Lactobacillales Enterococcaceae Enterococcus_E cecorum 0.020
bin_CaboVerde.16 Bacillota Bacilli Lactobacillales Lactobacillaceae Limosilactobacillus reuteri 0.000
bin_Denmark.56 Bacillota Bacilli Lactobacillales Lactobacillaceae Levilactobacillus brevis 0.000
bin_Malaysia.72 Bacillota Bacilli Lactobacillales Lactobacillaceae Limosilactobacillus reuteri 0.000
bin_CaboVerde.60 Bacillota Bacilli Lactobacillales Lactobacillaceae Ligilactobacillus agilis 0.001
bin_Malaysia.127 Bacillota Bacilli Lactobacillales Lactobacillaceae Ligilactobacillus agilis 0.001
bin_Denmark.61 Bacillota Bacilli Lactobacillales Lactobacillaceae Latilactobacillus sakei 0.003
bin_Malaysia.35 Bacillota Bacilli Lactobacillales Lactobacillaceae Ligilactobacillus animalis 0.003
bin_Malaysia.4 Bacillota Bacilli Lactobacillales Lactobacillaceae Ligilactobacillus ruminis 0.035
bin_CaboVerde.14 Bacillota Bacilli Lactobacillales Lactobacillaceae Lactobacillus acidophilus 0.067
bin_Aruba.18 Bacillota Bacilli Lactobacillales Streptococcaceae Streptococcus lutetiensis 0.001
bin_Brazil.22 Bacillota Bacilli Lactobacillales Streptococcaceae Streptococcus pasteurianus 0.001
bin_Denmark.117 Bacillota Bacilli Lactobacillales Streptococcaceae Lactococcus garvieae 0.001
bin_Denmark.113 Bacillota Bacilli Lactobacillales Streptococcaceae Streptococcus alactolyticus 0.002
bin_Brazil.69 Bacillota_A Clostridia Oscillospirales Acutalibacteraceae Clostridium_A leptum 0.097
bin_Denmark.93 Bacillota_A Clostridia Lachnospirales Anaerotignaceae Anaerotignum sp001304995 0.074
bin_CaboVerde.3 Bacillota_A Clostridia Clostridiales Clostridiaceae Clostridium_P perfringens 0.000
bin_Denmark.42 Bacillota_A Clostridia Clostridiales Clostridiaceae Clostridium_P perfringens 0.000
bin_Brazil.136 Bacillota_A Clostridia Clostridiales Clostridiaceae Clostridium thermobutyricum 0.041
bin_Brazil.3 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Faecalimonas sp900555395 0.026
bin_Brazil.89 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Faecalimonas sp900550235 0.034
bin_Denmark.63 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Anaerostipes hadrus 0.040
bin_Denmark.19 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Faecalimonas sp900556835 0.052
bin_Denmark.118 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Dorea_A longicatena 0.067
bin_Malaysia.108 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.081
bin_Brazil.113 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Anaerobutyricum sp900754855 0.095
bin_Denmark.50 Bacillota_A Clostridia Peptostreptococcales Peptostreptococcaceae Peptacetobacter sp900539645 0.053
bin_Brazil.107 Bacillota_A Clostridia Monoglobales_A UBA1381 CAG-41 sp900066215 0.016
bin_Brazil.159 Fusobacteriota Fusobacteriia Fusobacteriales Fusobacteriaceae Fusobacterium_B sp900541465 0.040
bin_Brazil.140 Fusobacteriota Fusobacteriia Fusobacteriales Fusobacteriaceae Fusobacterium_A sp900543175 0.046
bin_Malaysia.6 Fusobacteriota Fusobacteriia Fusobacteriales Fusobacteriaceae Fusobacterium_B sp900545035 0.088

8.3.3 Estimate - support plot

estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
    dplyr::select(genome,mean)

support <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
    dplyr::select(genome,support)

phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
  filter(genome %in% estimate$genome) %>%
  arrange(match(genome, estimate$genome)) %>%
  dplyr::select(phylum, colors) %>%
  unique() %>%
  arrange(phylum) %>%
  dplyr::select(colors) %>%
  pull()
Rows: 202 Columns: 2
── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (2): phylum, colors

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
inner_join(estimate,support,by=join_by(genome==genome)) %>%
    mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
    mutate(support=ifelse(mean<0,1-support,support)) %>%
    left_join(genome_metadata, by = join_by(genome == genome)) %>%
    mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
    ggplot(aes(x=mean,y=support,color=phylum))+
      geom_point(alpha=0.7, shape=16, size=3)+
      scale_color_manual(values = phylum_colors) +
      geom_vline(xintercept = 0) +
      xlim(c(-0.4,0.4)) +
      labs(x = "Beta regression coefficient", y = "Posterior probability") +
      theme_minimal()+
       theme(legend.position = "none")

8.4 Correlations

#Compute the residual correlation matrix
OmegaCor = computeAssociations(m)

# Refernece tree (for sorting genomes)
genome_tree_subset <- genome_tree %>%
        keep.tip(., tip=m$spNames) 


#Co-occurrence matrix at the animal level
supportLevel = 0.95
toPlot = ((OmegaCor[[1]]$support>supportLevel)
          + (OmegaCor[[1]]$support<(1-supportLevel))>0)*OmegaCor[[1]]$mean

matrix <- toPlot %>% 
      as.data.frame() %>%
      rownames_to_column(var="genome1") %>%
      pivot_longer(!genome1, names_to = "genome2", values_to = "cor") %>%
      mutate(genome1= factor(genome1, levels=genome_tree_subset$tip.label)) %>%
      mutate(genome2= factor(genome2, levels=genome_tree_subset$tip.label)) %>%
      ggplot(aes(x = genome1, y = genome2, fill = cor)) +
            geom_tile() + 
            scale_fill_gradient2(low = "#be3e2b",
                       mid = "#f4f4f4",
                       high = "#b2b530")+
            theme_void()

htree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************
vtree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************
#create composite figure
grid.arrange(grobs = list(matrix,vtree),
             layout_matrix = rbind(c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1)))

8.5 Predict responses

# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")

gradient = c("domestic","feral")
gradientlength = length(gradient)

#Treatment-specific gradient predictions
pred <- constructGradient(m, 
                      focalVariable = "origin", 
                      non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
            predict(m, Gradient = ., expected = TRUE) %>%
            do.call(rbind,.) %>%
            as.data.frame() %>%
            mutate(origin=rep(gradient,1000)) %>%
            pivot_longer(!origin,names_to = "genome", values_to = "value")
# weights:  9 (4 variable)
initial  value 101.072331 
final  value 91.392443 
converged

8.5.0.1 Element level

elements_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

community_elements <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(elements_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })

calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

element_predictions <- map_dfc(community_elements, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_elements[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

8.5.1 Positive associated functions at element level

# Positively associated

unique_funct_db<- GIFT_db[c(2,4,5,6)] %>% 
  distinct(Code_element, .keep_all = TRUE)


element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  filter(positive_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()
tinytable_6l6tn8kfxxyqtysuz2t3
trait mean p1 p9 positive_support negative_support Domain Function Element
B0103 0.008544525 0.0004896855 0.016616286 0.916 0.084 Biosynthesis Nucleic acid biosynthesis UDP/UTP
D0504 0.004425848 0.0003985114 0.009383968 0.938 0.062 Degradation Amino acid degradation Methionine
D0906 0.003599450 0.0001480281 0.007752579 0.926 0.074 Degradation Antibiotic degradation Fosfomycin
D0205 0.012498142 0.0019329087 0.023496427 0.944 0.056 Degradation Polysaccharide degradation Pectin
D0208 0.010205394 0.0017412487 0.019823177 0.938 0.062 Degradation Polysaccharide degradation Mixed-Linkage glucans
element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()
tinytable_4le7tbxfmlsou7saub0y
trait mean p1 p9 positive_support negative_support Domain Function Element
B0219 -0.003892460 -0.009445156 -1.598674e-04 0.047 0.953 Biosynthesis Amino acid biosynthesis GABA
B0214 -0.021086408 -0.039777480 -3.107959e-03 0.077 0.923 Biosynthesis Amino acid biosynthesis Glutamate
B0302 -0.004702329 -0.010402985 -4.605576e-04 0.036 0.964 Biosynthesis Amino acid derivative biosynthesis Betaine
B0310 -0.013024436 -0.024370029 -3.149799e-03 0.048 0.952 Biosynthesis Amino acid derivative biosynthesis Tryptamine
B0303 -0.011613334 -0.021373921 -2.413267e-03 0.064 0.936 Biosynthesis Amino acid derivative biosynthesis Ectoine
B0804 -0.016684242 -0.031511256 -3.389357e-03 0.052 0.948 Biosynthesis Aromatic compound biosynthesis Dipicolinate
B0603 -0.016440795 -0.031679497 -2.316252e-03 0.070 0.930 Biosynthesis Organic anion biosynthesis Citrate
B0601 -0.009340160 -0.018754997 -7.868225e-04 0.084 0.916 Biosynthesis Organic anion biosynthesis Succinate
B0401 -0.012342178 -0.023948940 -4.468990e-04 0.094 0.906 Biosynthesis SCFA biosynthesis Acetate
B0709 -0.002061541 -0.003660689 -5.929844e-04 0.040 0.960 Biosynthesis Vitamin biosynthesis Tocopherol/tocotorienol (E)
D0517 -0.004458863 -0.007873548 -1.198957e-03 0.030 0.970 Degradation Amino acid degradation Ornithine
D0508 -0.003160979 -0.007357970 -1.369150e-04 0.088 0.912 Degradation Amino acid degradation Lysine
D0903 -0.003870200 -0.009394756 -1.560699e-04 0.047 0.953 Degradation Antibiotic degradation Cephalosporin
D0908 -0.015263907 -0.027752320 -3.008654e-03 0.071 0.929 Degradation Antibiotic degradation Macrolide
D0601 -0.008936709 -0.017382607 -1.878995e-03 0.042 0.958 Degradation Nitrogen compound degradation Nitrate
D0611 -0.003870200 -0.009394756 -1.560699e-04 0.047 0.953 Degradation Nitrogen compound degradation Phenylethylamine
D0603 -0.001894066 -0.003700352 -3.047402e-04 0.053 0.947 Degradation Nitrogen compound degradation Urate
D0610 -0.002980443 -0.004927790 -8.438420e-04 0.059 0.941 Degradation Nitrogen compound degradation Methylamine
D0606 -0.005546897 -0.011162503 -6.350907e-04 0.073 0.927 Degradation Nitrogen compound degradation Allantoin
D0612 -0.001551679 -0.002998413 -8.337605e-05 0.089 0.911 Degradation Nitrogen compound degradation Hypotaurine
D0801 -0.001834932 -0.002282080 -9.835430e-05 0.009 0.991 Degradation Xenobiotic degradation Toluene
D0802 -0.001834932 -0.002282080 -9.835430e-05 0.009 0.991 Degradation Xenobiotic degradation Xylene
D0807 -0.003989919 -0.008479592 -6.243554e-04 0.057 0.943 Degradation Xenobiotic degradation Catechol
D0817 -0.004699653 -0.010697649 -5.137227e-04 0.060 0.940 Degradation Xenobiotic degradation Trans-cinnamate

8.5.2 Negatively associated functions at element level

element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()
tinytable_4v8q9k8cmxg739les76w
trait mean p1 p9 positive_support negative_support Domain Function Element
B0219 -0.003892460 -0.009445156 -1.598674e-04 0.047 0.953 Biosynthesis Amino acid biosynthesis GABA
B0214 -0.021086408 -0.039777480 -3.107959e-03 0.077 0.923 Biosynthesis Amino acid biosynthesis Glutamate
B0302 -0.004702329 -0.010402985 -4.605576e-04 0.036 0.964 Biosynthesis Amino acid derivative biosynthesis Betaine
B0310 -0.013024436 -0.024370029 -3.149799e-03 0.048 0.952 Biosynthesis Amino acid derivative biosynthesis Tryptamine
B0303 -0.011613334 -0.021373921 -2.413267e-03 0.064 0.936 Biosynthesis Amino acid derivative biosynthesis Ectoine
B0804 -0.016684242 -0.031511256 -3.389357e-03 0.052 0.948 Biosynthesis Aromatic compound biosynthesis Dipicolinate
B0603 -0.016440795 -0.031679497 -2.316252e-03 0.070 0.930 Biosynthesis Organic anion biosynthesis Citrate
B0601 -0.009340160 -0.018754997 -7.868225e-04 0.084 0.916 Biosynthesis Organic anion biosynthesis Succinate
B0401 -0.012342178 -0.023948940 -4.468990e-04 0.094 0.906 Biosynthesis SCFA biosynthesis Acetate
B0709 -0.002061541 -0.003660689 -5.929844e-04 0.040 0.960 Biosynthesis Vitamin biosynthesis Tocopherol/tocotorienol (E)
D0517 -0.004458863 -0.007873548 -1.198957e-03 0.030 0.970 Degradation Amino acid degradation Ornithine
D0508 -0.003160979 -0.007357970 -1.369150e-04 0.088 0.912 Degradation Amino acid degradation Lysine
D0903 -0.003870200 -0.009394756 -1.560699e-04 0.047 0.953 Degradation Antibiotic degradation Cephalosporin
D0908 -0.015263907 -0.027752320 -3.008654e-03 0.071 0.929 Degradation Antibiotic degradation Macrolide
D0601 -0.008936709 -0.017382607 -1.878995e-03 0.042 0.958 Degradation Nitrogen compound degradation Nitrate
D0611 -0.003870200 -0.009394756 -1.560699e-04 0.047 0.953 Degradation Nitrogen compound degradation Phenylethylamine
D0603 -0.001894066 -0.003700352 -3.047402e-04 0.053 0.947 Degradation Nitrogen compound degradation Urate
D0610 -0.002980443 -0.004927790 -8.438420e-04 0.059 0.941 Degradation Nitrogen compound degradation Methylamine
D0606 -0.005546897 -0.011162503 -6.350907e-04 0.073 0.927 Degradation Nitrogen compound degradation Allantoin
D0612 -0.001551679 -0.002998413 -8.337605e-05 0.089 0.911 Degradation Nitrogen compound degradation Hypotaurine
D0801 -0.001834932 -0.002282080 -9.835430e-05 0.009 0.991 Degradation Xenobiotic degradation Toluene
D0802 -0.001834932 -0.002282080 -9.835430e-05 0.009 0.991 Degradation Xenobiotic degradation Xylene
D0807 -0.003989919 -0.008479592 -6.243554e-04 0.057 0.943 Degradation Xenobiotic degradation Catechol
D0817 -0.004699653 -0.010697649 -5.137227e-04 0.060 0.940 Degradation Xenobiotic degradation Trans-cinnamate
positive <- element_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- element_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

table <- bind_rows(positive,negative) %>%
  left_join(unique_funct_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) 

table %>%
  mutate(Element=factor(Element,levels=c(table$Element))) %>%
  ggplot(aes(x=mean, y=fct_rev(Element), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.5.2.1 Function level

functions_table <- elements_table %>%
    to.functions(., GIFT_db) %>%
    as.data.frame()

community_functions <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(functions_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })
#max-min option
calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

function_predictions <- map_dfc(community_functions, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_functions[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

# Positively associated
function_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  tt()
tinytable_2ivmnqw185r7pbswxgd3
trait mean p1 p9 positive_support negative_support
D02 8.144008e-03 -0.003585743 0.0201994472 0.820 0.180
B08 7.729591e-03 -0.003619195 0.0188190890 0.782 0.218
B01 1.079193e-03 -0.006174439 0.0083360499 0.643 0.357
S01 4.929735e-04 -0.012595581 0.0138898138 0.556 0.444
D01 4.722996e-05 -0.005072307 0.0047107372 0.430 0.570
B09 8.550675e-05 -0.000537705 0.0005261702 0.375 0.625
# Negatively associated
function_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  tt()
tinytable_3dgk3vvu4nkgtl5u74ji
trait mean p1 p9 positive_support negative_support
D08 -1.128130e-03 -0.0021383271 -0.0001797254 0.046 0.954
B03 -1.050395e-02 -0.0186128263 -0.0028600447 0.056 0.944
D06 -3.041585e-03 -0.0069243804 0.0001474005 0.112 0.888
B04 -8.359605e-03 -0.0175260241 0.0009514529 0.133 0.867
B06 -6.693422e-03 -0.0163287939 0.0032004529 0.184 0.816
D07 -1.184844e-02 -0.0305629696 0.0040330573 0.194 0.806
D05 -1.464921e-03 -0.0078777280 0.0038694072 0.214 0.786
D03 -3.992202e-03 -0.0126804243 0.0032883113 0.224 0.776
S03 -9.018483e-03 -0.0330709641 0.0152032371 0.250 0.750
B02 -3.531226e-03 -0.0122475416 0.0047519621 0.270 0.730
D09 -1.973050e-03 -0.0085424585 0.0046762375 0.307 0.693
S02 -4.605229e-03 -0.0136699132 0.0030797554 0.323 0.677
B07 -3.729566e-03 -0.0158138065 0.0091826160 0.330 0.670
B10 -8.894201e-06 -0.0003238718 0.0002411792 0.484 0.516
positive <- function_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- function_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.02,0.02)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")